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Light from 'point sources' such as supernovae is observed with a beam width of order of the 
sources' size - typically less than 1 AU. Such a beam probes matter and curvature distributions 
that are very different from coarse-grained representations in N-body simulations or perturbation 
theory, which are smoothed on scales much larger than 1 AU. The beam typically travels through 
unclustered dark matter and hydrogen with a mean density much less than the cosmic mean, and 
through dark matter halos and hydrogen clouds. Using N-body simulations, as well as a Press- 
Schechter approach, we quantify the density probability distribution as a function of beam width 
and show that, even for Gpc-length beams of 500 kpc diameter, most lines of sight are significantly 
under-dense. From this we argue that modelling the probability distribution for AU-diameter beams 
is absolutely critical. Standard analyses predict a huge variance for such tiny beam sizes, and 
nonlinear corrections appear to be non-trivial. It is not even clear whether under-dense regions lead 
to dimming or brightening of sources, owing to the uncertainty in modelling the expansion rate which 
we show is the dominant contribution. By considering different reasonable approximations which 
yield very different cosmologies we argue that modelling ultra-narrow beams accurately remains 
a critical problem for precision cosmology. This could appear as a discordance between angular 
diameter and luminosity distances when comparing SN observations to BAO or CMB distances. 



I. INTRODUCTION: ON NARROW BEAMS 

Supernovae la (SNIa) observations play a critical role 
in the evidence for a nonzero cosmological constant (see 
[T] and references therein). SNIa are effective standard 
candles (we think their intrinsic luminosity can be cali- 
brated from their light curve). In the standard cosmo- 
logical model, their observed luminosity is used to in- 
fer their luminosity distance (or equivalently magnitude) 
by assuming that the geometry of the universe is well- 
described by a Friedmann-Lemaitre (FL) background ge- 
ometry that describes the universe smoothed on large 
scales. Of course, photons do not propagate in the ge- 
ometry of a smooth FL spacetime but in the lumpy uni- 
verse: a beam mostly propagates in underdense regions 
between clustered matter (overdense islands of matter). 
The problem of quantifying the effects of propagation in 
an inhomogeneous universe was first addressed, as far as 
we are aware, independently by Zel'dovich |2] and Feyn- 
man [3]. Zel'dovich introduced the empty-beam approxi- 
mation to deal with light rays propagating in vacuum and 
this was extended to the case of a partially-filled beam by 
[4j. These results later came to be known as the Dyer- 
Roeder approximation - see below. Zel'dovich's insight 
also led to other work on the problem in the 1960s [SHE] ■ 

Light propagation in inhomogeneous spacetimes gives 



rise both to distortion of images and magnification of 
some images, because of gravitational lensing; in com- 
pensation, most images are demagnified. These effects 
induce, in particular, a dispersion of the observed SNIa 
luminosities and hence an extra scatter in the Hubble 
diagram [5HT5]. "Precision cosmology" within the stan- 
dard approach could be compromised by the effects of 
lensing on the interpretation of SNIa data - and thus it 
is crucial to characterise the magnitude of these effects 
precisely. A related key question is "what physical and 
angular sizes are relevant in estimating these effects on 
SNIa observations?" 

A perturbative approach (i.e. with light propagating 
in a perturbed FL spacetime) shows that the dispersion 
due to large-scale structure becomes comparable to the 
intrinsic dispersion for redshifts z > 1 [14j . However, 
since the matter fluctuations responsible for the magnifi- 
cation of the SNIa also induce a shearing of the images of 
background galaxies, this dispersion can actually be cor- 
rected [m [16] . A similar idea was pursued in the context 
of gravitational wave sirens [T7]. Nevertheless, a consid- 
erable fraction of the lensing dispersion arises from sub- 
arcminute scales, which are not probed by shear maps 
smoothed on arcminutc scales [18]. 

To estimate the dispersion induced by inhomogeneities, 
one first needs to determine the typical size of the 
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geodesic bundle associated with SNIa. The typical ob- 
servational aperture is of order 1", whereas the beam is 
actually much thinner: 0(1) AU for a source at redshift 
z ~ 1, i.e. an aperture of /? ^ 10^^". This is typically 
smaller than the mean distance between any massive ob- 
jects (galaxies, stars, H clouds, small dark matter halos) 
- and on a scale where the fluid continuum model may 
not be suitable any more. Thus the beam propagates 
in preferentially low density regions with rare encounters 
of gravitationally collapsed, high density patches (halos) 
resulting in highly inhomogeneous geometry. 

By contrast, the standard approach implicitly uses a 
perturbative analysis and a fluid continuum model by 
treating the beam as propagating in the background and 
perturbed FL geometry. From this viewpoint, it is sur- 
prising that the standard analysis of SNIa data leads to 
a consistent result, in particular with other cosmological 
probes. Is a smooth cosmological model a good descrip- 
tion of our universe, and in particular for interpreting 
data such as SNIa? If so, can we understand clearly why 
this is the case? 

Standard perturbation theory reveals there are prob- 
lems. If the angular diameter distance as a function of 
redshift in a perturbed FL model differs from the back- 
ground by k(/3) where /? is the angular scale of observa- 
tion (see below for definitions), then the variance of k, 
scales as, 

{^m''^^lQ-^a,nr'zr{ff''^'\ (1) 

as derived in [121 120], using linear perturbation theory 
(n is the spectral index of the power spectrum of den- 
sity fluctuations and Zg the redshift of the source). This 
estimate was confirmed in |21| considering the nonlinear 
evolution of the power spectrum, but on scales below 10 
arcmin, (k^(/3)) increases more steeply than the theoret- 
ical expectation of the linear theory and is 2 to 3 times 
higher. Note also that IT]) implies that (k^(/3))^/^ be- 
comes of order 2 x IQ^ on an angular scale of order 
1 arcmin, so that the variance becomes much larger than 
unity on the typical angular size of the ray bundle for 
SNIa. So large-scale structure induces a stochastic dis- 
persion of the luminosity distance which is difficult to 
quantify with standard techniques for narrow beams. As 
soon as one goes down to much smaller scales, inhomo- 
geneities also induce a systematic shift away from the 
background, since one cannot neglect the higher order 
terms. This is much more serious. The magnification of 
a source behaves as ^ ~ 1-f 2K-f 3K^-f |7p-|-. . ., where 7 is 
the shear of the source (defined below), and so the mean 
of the magnification is {^) ^ 1 + {3k^ + I7P) + ... 7^ 1. 
Thus, if the variance is large on small angular scales, we 
expect the overall shift to be significant - potentially of 
order unity or larger - on these scales too. Modelling this 
shift accurately is critical for interpreting SNIa observa- 
tions correctly. 

To estimate the effect of the inhomogeneities on 
smaller scales, one needs to provide a better description 



of the matter distribution. Attempts to include a uniform 
component, high density halos and low density zones (fil- 
aments and voids) have been proposed [311 US] j but none 
go down to the required scale. 

Narrow light bundles travel large distances (> 
100 h^^Mpc) with a very low probability of encountering 
dark matter halos of substantial mass, which we quantify 
in the following section. The cuspy density profiles of the 
halos additionally reduce the probability of a bundle to 
cross the central, high density regions. Thus the bun- 
dles are subject mainly to Weyl focussing (i.e. induced 
by the gradient of the gravitational potential). These 
lightrays are expected to be demagnified compared to 
lightrays propagating in a FL spacetime of mean mat- 
ter density. When one averages such ray bundles over 
the whole sky, this dimming is compensated by a small 
number of ray bundles that encounter very large density 
inhomogeneities , which then results in high focussing for 
those directions and a magnification. The general aver- 
aging argument was put forward by Weinberg [24j . 

A similar argument can be employed to determine the 
average density encountered by a light bundle from a su- 
pernova. The probability of a light bundle encountering 
massive halos decreases with halo mass. Massive galax- 
ies, groups and clusters (> 10^^ h~^MQ) are rarely en- 
countered, yet, they comprise ~ 50% of total mass within 
the universe. Thus, if we measure SNIa in typical di- 
rections in the sky, we observe them in directions where 
the density of matter encountered by the relevant ray 
bundles may be expected to be less than the cosmologi- 
cal average (for almost all directions are of this nature) . 
This argument does not apply to the much larger angular 
scales relevant to measurements of the BAO and CMB 
peaks. These beams do encounter sufficient matter to on 
average correspond to the overall cosmological density, as 
argued by [23]. In these circumstances, one expects not 
only an extra dispersion in the Hubble diagram, induced 
by the spatial inhomogeneity of the intervening medium, 
but also a systematic shift, induced by an observational 
selection effect, which may well be significant. 

The goal of this article is to investigate these effects. 
We first discuss modelling the matter distribution in the 
real universe. Then we consider the general relativistic 
problem of light rays in an arbitrary spacetime, with a 
focus on how a light beam reacts to inhomogeneities com- 
pared to a smooth spacetime. We then consider some 
different approximations used to model narrow beams, 
such as perturbed FL models, the Dyer-Roeder approx- 
imation, as well as presenting two new approximations. 
Finally, we consider the problem of how Weyl focussing 
by many point sources is converted into Ricci focussing 
associated with a smooth matter distribution. 



II. THE MATTER DISTRIBUTION 

According to the most accepted variants of the ACDM 
paradigm gravitationally collapsed structures span a 
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FIG. 1: The top panel shows the total fraction of cosmic 
mass which is locked in halos above a given mass. The 
Press-Schechter model predicts 50% of the total mass is 
locked in halos above 3 x 10^^ h~^MQ. The second panel 
gives the number of halos per mass bin per /i~^Mpc'^. 
The third panel displays the average number of halos per 
mass bin encountered by an infinitesimal thin light ray 
per unit length. The bottom panel indicates the proba- 
bility that a ray intersects with at least one halo within 
a given mass bin on a distance of 1 /i~^Mpc, which is 
equivalent to 1 — Ph-cc, where Pfrco is the probability of 
not hitting a halo within that mass bin on a length of 
1 h-^Mpc. 



mass range from 10~* Mq (determined by the free- 
streaming scale of the CDM particle) to 10^^ Mq. In 
addition to the matter bound in collapsed structures a 
smooth component is expected which is not bound to 
any structure. The question we are interested in here 
is, how these structures affect the propagation of light 
bundles travelling through this inhomogeneous universe. 

As first attempt to answer to this question we employ 
an analytic approach based on a Press-Schechter (PS) 
model [25]. Press and Schechter derived an analytical 
expression for the cumulative mass function, i^(Mhaio), 
which gives the fraction of mass locked in halos above 
a given mass. Integrating their formula over the whole 



mass spectrum yields one. Thus the PS model predicts 
that all matter is bound in gravitationally collapsed ha- 
los, with masses ranging from zero to infinity. With the 
increasing dynamical range of N-body simulations, quan- 
titative differences with this model have become apparent 
[26H33j. However, the difference between analytical and 
N-body predictions depends on the definition of a halo in 
the simulations. Recently, 'M!,'351 showed that using 'dy- 
namical masses' for N-body halos yields good agreement 
at least for low redshifts. 

The top panel of Fig. [T] shows our computation of 
-F(Mhaio) based on the 5- Year WMAP data [36]. The 
derivative of F determines the number density of halos, 
-^haio, as a function of mass (second panel). According to 
spherical collapse theory the radius of a halo with mass, 
Afhaio, is given by: 

/ SMhalo /„^ 

''halo = ;p-r (2) 

where pcrit is the critical density of the universe and Ac = 
95 (for cosmological parameters chosen above). 

The average number of halos, AThaio, with mass A/haio 
encountered by a single infinitesimal ray per unit length 
is the product of A'haio and the surface area of the halo. 

-'^halo = A^haloTrr^a^jQ (3) 

^haio is shown in the third panel from the top. The 
average number of halos with masses about 10^^ h~^MQ 
encountered by a light ray is ^ 0.001. Thus, on average a 
single light ray of light hits one 10^^ H^^Mq halo every 
1000 h-^Mpc. 

The bottom panel displays the probability that a ray 
intersects with at least one halo of a given mass while 
travelling a distance of 1 /i~^Mpc, which is equivalent 
to 1 — Pfree, whcrc Pfree is the probability of passing 
1 h~^Mpc freely, i.e., without encountering a single halo 
of that mass. To compute Pfree we subdivide the volume 

— 1/3 

into cubes with I/A^id — A'j^aio '^^ ^ ^^^^ ^^^'^ assume 
that halos within a given mass bin are distributed quasi 
homogeneously, i.e, only one single halo is placed within 
each cube. The position of the halo within the cube is 
chosen randomly. With these assumptions Pfrcc can be 
approximated by 

Pfr.. ^ {1 - nrUN^^f''' . (4) 

We conclude that a light ray travelling 1 h~^Mpc 
through the present, nonlinear, cosmic density field en- 
counters with almost 100% certainty several halos with 
masses below 10~^ h~^MQ. Since the PS model does 
not include a smooth component the predictions for this 
mass range must be corrected accordingly. On the other 
side, the probability to hit at least one 10"'^^ hr^MQ halo 
is of the order of 0.1%. These results suggest we must 
discuss the effects of small and large scale structure on 
the light propagation separately. 
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A. Effects of small scale structure 

Lensing can discriminate between a diffuse and smooth 
component (a gas of microscopic particles) and one 
of macroscopic massive objects (gravitationally bound), 
and has been used [371 131] to probe the nature of dark 
matter on galactic scales. The two components can be 
characterized by a mass scale, defined by the fact that 
the projected density be smooth on a scale of order the 
angular size of the source. This gives [37] 



M, 



2 X IQ-'^^Mp.h^ 



As 
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(5) 



where As is the physical size of the source. 

Another important mass scale is set by the requirement 
that the angular size of the source, /3 — Xs/ Da{zs), is 
smaller than the Einstein angular radius 0e so that it 
can be considered as a true point source [37) : 
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If M < M^, the component can be considered as diffuse 
on the scale of the ray bundle. If M* < M < Mpomt there 
will be very few high magnification events with most of 
the lines of sight being demagnified, according to the 
standard lensing paradigm (see below). These two com- 
ponents affect the probability distribution function for 
the magnification. In the extreme case where the matter 
is composed only of macroscopic point-like objects, then 
most high-redshift SNIa would appear less bright than in 
a universe with the same density distributed smoothly, 
with some very rare events of magnified SNIa [37 [ l40 l |4T] . 
It has been argued [35] that the magnification of high- 
redshift SNIa can be a powerful discriminator of the na- 
ture of dark matter. In particular, based on numerical 
simulations, a few hundred SNIa at z ~ 1 could allow a 
20% determination of the fraction of matter in compact 
objects. 

The dispersion of SNIa data due to lensing has been 
estimated in various ways, but a complete analysis may 
require us to go down to scales where our knowledge of 
the distribution of matter is very poor. In particular, it 
is important to know the amount of diffuse matter com- 
pared to the amount of matter in small compact halos. 
Knowledge of the spatial distribution of the halos is re- 
quired to determine the dispersion of the observations, 
keeping in mind that one also expects a bias since most 
lines of sigh are demagnified. 

The minimum mass of gravitationally bound structures 
is determined by the nature of the CDM particle itself. 
Its mass induces a free streaming scale which in turn gives 
rise to a mass threshold below which no gravitationally 
bound structure can form (in contrast to the original PS 
approach where there is no such cut off mass). Cur- 
rently, neutralinos are the most promising CDM candi- 
dates. The lightest neutralino, with m ^ 100 GeV, is 



favoured, as it is both weakly interacting and stable (e.g. 
[43]). Its free-streaming scale is ~ Q.7pc, with a corre- 
sponding minimum halo mass Mfg ~ 1O~*M0. 

In principle, cosmological N-body simulations can be 
employed to determine the total amount of mass in the 
smooth and the halo component. Knowing the funda- 
mental properties of the CDM particle, the initial condi- 
tions can be derived and propagated to the current epoch. 
However, the dynamical range required for this approach 
exceeds current computational resources. 

One strategy to circumvent the computational limita- 
tions is to enclose a small region of very high resolution 
within a larger but lower resolution simulation . Using 
this technique [H] found that at z = 26 the mass func- 
tion is steep, dn{M)/dM cx down to Mfs. At that 
time about 5% of the mass in the high resolution region 
has collapsed into gravitationally bound halos (see Fig. 
3 in [331). Due to technical limitations this simulation 
has not been run further than z « 26. 

Another strategy to determine the total mass locked in 
halos is via excursion set theory [3S], which propagates 
density perturbations stochastically to generate the halo 
mass functions. For a ACDM model with 100 GeV neu- 
tralinos, 75-80% of matter is locked in halos at z < 1 [46) . 
The remaining 20-25% is smoothly distributed without 
being associated with any collapsed structure. 



B. Effects of large scale structure 

Current N-body simulations provide a very reliable 
picture of the cosmic large scale structure. However 
small scale structure can only be resolved down to the 
given mass resolution limit of these simulations (currently 
~ 10^ -^o)- Analytical approaches, like those based on 
PS models, are not affected by mass resolution issues but 
are much more sketchy by nature. In the following we will 
investigate both approaches. 



1. Analytical approach 

Based on the PS model discussed above one can de- 
termine the probability distribution (PDF) of densi- 
ties averaged along infinitesimally thin light rays. For 
that purpose we model the mass distribution of halos 
by four bins with average masses from 10^ h~^MQ to 
3 X 10^* h~^M0. The mass bins are chosen in such a way 
that the one dimensional number density, A^id, increases 
by a factor of 10 for each subsequent (decreasing) mass 
bin. The remaining mass contained within halos below 
the smallest mass bin is assumed to be homogeneously 
distributed. According to the PS model 15% of the total 
gravitational matter is found in halos below 10^ /i^^M©. 

The probability of a ray encountering a given total 
number of halos is computed as product of the various 
binomial coefficients and probabilities for hitting the par- 
tial number of halos per mass bin. 
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FIG. 2: Probability distribution of the averaged densities 
encountered by a single light ray of infinitesimal width for 
different path lengths based on a Press-Schechter model. 

This approach allows us to compute the PDFs as a 
function of the path length as shown in Fig. [2] The dis- 
tribution for a path length of 100 /i~^Mpc peaks at 0.4 
times the mean density, Pmcan- For longer path lengths 
the peak shifts towards Pmcan- But even for a path length 
of 1000 /i^^Mpc (corresponding to a redshift of z ~ 0.25) 
the distribution peaks significantly below Pmean- It is 
worth noting that by construction the integral of p dP 
from zero to infinity equals unity. A peak below Pmoan 
requires a high density tail for counterbalance. 

The model presented here does not include halo clus- 
tering (inherent to such kind of approaches) and assumes 
the mass contained in objects below 10^ h~^MQ to be 
distributed smoothly (to reduce computational cost). 
These shortcomings let us abstain from a quantitative 
interpretation at this point. But we can clearly see that 
the majority of light rays encounters averaged densities 
below the mean and that the shape of the PDFs is a 
function of path length. 

2. Numerical approach 

Numerical simulations provide detailed insight into the 
nonlinear matter evolution inaccessible to analytical de- 
scriptions. However, they are inevitably limited in their 
mass resolution. The resolution is proportional to the 
simulation volume and inversely proportional the num- 
ber of phase space elements (particles) used. For cos- 
mological applications generally large volumes are desir- 
able. The number of particles is limited by the available 
computational power. Particles of current state-of-the- 
art cosmological simulations have masses of a few times 
10^ Mq. For comparison, the total mass, assuming ho- 
mogeneous density distribution, contained within the vol- 
ume of a light beam of 1 AU diameter and 1000 Mpc 



length (corresponding to z « 0.25) is '-^ 10 ^ Mq. 

In this section we compute the PDFs of the mean den- 
sities within long and narrow 'beams' based on a set of 
publicly available N-body simulations, namely the Bol- 
shoi EZj, the Millennium [iS' and the MultiDark Rl ^ 
Simulation. (Descriptions of the data bases are given 
in [SOI EI]). These simulations compute the CDM dis- 
tribution within cubes of 250^, 500^ and lOOO^ h'^Mpc^ 
with mass resolution of 1.3 x 10^, 8.6 x 10^ and 8.7 x 
10^ h~^MQ, respectively. These mass resolutions only 
allow the determination of the densities within beams 
wider than a few tens of kpc which is many orders of 
magnitude larger than the expected diameter of a light 
beam from a supernova. Nevertheless the results derived 
here can give basic insight into the mean densities within 
the volume of the light beams from distant SNla. 

The right panel of Fig. |3] shows the PDFs of the mean 
densities within beams of 500 h^^kpc diameter as a func- 
tion of their length. The dashed lines based on the Multi- 
Dark simulation are shown as a consistency check. They 
are expected to coincide with the Bolshoi results but are 
more affected by Poisson noise due to the lower mass res- 
olution. The shape of the distributions is similar to those 
based on the PS model (Fig. [2|. Independent of length 
all distributions peak below the mean density. The cu- 
mulative probability for a mean density below the cosmic 
mean for the 100, 250, 500 and 1000 h-^Mpc beams is 
75%, 71%, 68% and 65%, respectively The large proba- 
bility of an averaged density below cosmic mean is coun- 
terbalanced by unlikely high density encounters. Investi- 
gations based on the Millennium simulation confirm these 
results. The middle panel of Fig. [3] shows the PDFs for 
different diameters. There may be an indication that 
with decreasing diameters the location of the peak con- 
verges towards values close to 0.5 x pmcan, but we are un- 
able to probe beam sizes below ^ 50 /i~^kpc. The shape 
of the distribution, and in particular that the location 
of the peak is below unity, is preserved independent of 
diameter and is expected to hold for much smaller diam- 
eters. The right hand panel of Fig. |3] displays the PDFs 
based on cubic volumes, i.e. the densities are measured 
in cubes rather than beams (left panel) while the vol- 
ume remains the same. The overall shape of the PDFs is 
very similar to those shown on the left but the location 
of the peak is shifted towards smaller densities. Obvi- 
ously, the geometry of the 'test volumes' has an impact 
on the PDFs. The PDFs can not be accurately deter- 
mined without incorporating their spatial information of 
the large scale density distribution. 

The picture arising from the N-body simulations is con- 
sistent with the PS model. The probability that a light 
beam from a supernova encounters an average density 
less then cosmic mean is larger than 50%. The exact 
value depends on light path length. With shorter dis- 
tances the cumulative probability of sampling densities 
below cosmic mean increases. This effect may be suffi- 
cient to induce biases in the luminosity distance relation. 
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FIG. 3: Left panel: Probability distribution of averaged densities within long narrow (0.5 h~^Mpc diameter) beams 
for different path lengths using different simulations (indicated). For short path lengths (e.g., ~ 100 /i^^Mpc) most 
beams encounter low densities which is counterbalanced by comparatively few beams which encounter much higher 
densities. With increasing path length the peak of the PDF approaches the cosmic mean, but even for a beam 
length of 1 h~^Gpc the peak occurs significantly below the mean density. Independent of length the average density 
encountered by a sufficiently large number of beams is equal to the cosmic mean density, i.e., the density weighted 
integrals for all curves shown above yields cosmic mean density. Middle panel: PDF for beams of the same length 
(250 h~^Mpc) but different diameter (indicated). As the beam becomes narrower the PDF broadens and the location 
of the peak tends to shift to slightly smaller densities. Right panel: Same as left panel except here the density is 
measured within cubes which have the same volume as the long and narrow beams. Note the striking difference 
between the PDFs: the peak position and widths change when comparing tubes to cubes, and the power law tail 
which is prominent for beams is no longer present. 



III. LIGHT PROPAGATION 



(u^'u^t — —1), the rcdshift is defined by 



From a theoretical point of view, the effects of matter 
inhomogeneities can be described by the geodesic devia- 
tion equation, which describes the evolution of a bundle 
of geodesies x^{v, s), where v is the affine parameter and 
s labels the geodesies. The past lightcones of the central 
observer are given by w = const, where w is the phase. 
Then /c^ ~ d^w, so these curves are irrotational null 
geodesies: 



k^k^ = 0, fc^V.fc^ = 0, V[^fc,] = 0. 



(7) 



The connecting vector 77'^ = dx'^/ds relates neigh- 
bouring geodesies with tangent vector fc'^ = da;'^/dw to 
an arbitrary reference geodesic of the bundle, x'^{v) — 
x^{v,0), giving the distance between neighbouring 
geodesies and hence the physical size and shape of the 
bundle as one follows it down into the past. The con- 
necting vector can always be chosen such that fc^r;^ = 
and it evolves according to the geodesic deviation equa- 
tion: 



k^k'^VaVpf]" = R''^c.pk''k"r]'^ 



(8) 



l + z{v) 



{k^U'')y 



where the past-directed photon four-momentum is 



fc''==(l + z)(-ii^ + e^), e^u^ = 0, e'^ep 



1. 



(9) 



(10) 



Here is the spatial direction of observation, and the 
spatial direction of propagation is 71'' = — e^. The affine 
parameter increases monotonically along each ray and co- 
incides in an infinitesimal neighborhood of the observa- 
tion point with the Euclidean distance in the rest frame 
of u^{0). Note that while it depends on the 4- velocity 
u^ {0) of the observer, it does not depend on the 4- velocity 
u^ {x°'{s)) of the observed source. 

The screen space at each point along a ray is in the 
observer's rest space and orthogonal to the ray direc- 
tion. It is spanned by unit vectors (a = 1,2), with 
g^uTi'^n'^ = Sab and n^w^ = ni^k^ = 0, that are parallel 
transported along the ray (/c^V^n^ = 0). We can choose 
the connecting vector to lie in the screen space, so that^ 



This equation describes the change of shape of the bun- 
dle. 

For fundamental observers with four-velocity 



^ This is the Sachs basis, unique up to transformations 
r[if,(a)nj^ +Paki^, where rabi^t) is a rotation through angle o, 
and Pa are constants. 
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d2 



Va = 'R-abV , 



(11) 



where TZab = Rp.vafik'^k^n'^n^ is the screen projection of 
the Riemann tensor. We write 



T^ab — 



$00 
$00 



-RcvJ/q Im^'o 
Im ^0 Re ^0 



(12) 



with 



$n 



*o = -lc^,o.pm''k''m^kP, (13) 



where m'^ = Tn\ — m!^. The Einstein equations give 
R^^kt^k" = STrGT^i.Ai^fc'', where T^,. is the total energy- 
momentum tensor, 



(p + p)u^u^ + pgi_,y + TT^^ + q^u^ + q^i 



(14) 



Here 7r^,y is the anisotropic stress and (7^ is the momen- 
tum density. (For a perfect fluid tt^^ = = g^; for 
more general fluids, we can always choose = 0, corre- 
sponding to the frame where comoving observers see no 
momentum flux). Then we find 

$oo(v) = -4ttG[1+z{v)]^ {p + P + 2q^e^ + Tr^.e^e'') 

(15) 

Note that a cosmological constant A makes no contribu- 
tion to $00 • 



The linearity of (111 implies that 

dry'' 



dv 



(16) 



v=0 



where is the Jacobi map. By ( 11 1, we have the Jacobi 
matrix equation 

du^ dv 

(17) 

This second-order linear equation can be rewritten as a 
first-order nonlinear equation: 



dv 



by defining the deformation matrix 
d 



dv 



(18) 



(19) 



The Jacobi map V^b or equivalently the deformation ma- 
trix iS"b are the central quantities to describe the distor- 
tion of the geodesic bundle. The deformation matrix is 
usually decomposed as^ 







d-2 \ 






-<^1 1 



(20) 



Recall that the null rotation Vj^/c^j vanishes since fc^ = 9,jUi. 



which defines the optical scalars 9 (null expansion) and 
a = ai + ii72 (null shear). These satisfy the Sachs equa- 
tions [S3] 



d^ 

dv 



do- 

dv 



29a 



$00, 

'J'o, 

1 



Vak^Vk'' - 9' 



(21) 
(22) 
(23) 



The evolution of a ray bundle can then be discussed in 
terms of Ricci focussing ($00) and Weyl focussing (vl/o)- 
The first is generated by matter inside the beam [see 
(15)] while the second derives from matter outside the 
beam, which can generate a non-vanishing Weyl tensor 
inside the beam. This distinction leads to the problem 
raised by Zel'dovich [2] and Feynman |3], and posed in 
terms of the curvature tensor by Bertotti fS] : if the mat- 
ter of the universe is clustered in massive galaxies, the 
bundle propagates almost exclusively in vacuum, or at 
least in underdense regions, and is thus mostly subject 
only to the Weyl focussing; by contrast, the cosmological 
effect is modelled using a homogeneous fluid which gen- 
erates only Ricci focussing (the Weyl tensor vanishes in 
FL spacetime). Dyer and Roeder [SJ (see also [55H57| ) 
effectively reproduced Zel'dovich's idea and proposed an 
ansatz to model the propagation in regions with no inter- 
galactic medium. Weinberg [24] disputed this model, ar- 
guing that multiple Weyl deflections by individual masses 
average to mimic the Ricci effect of a fluid with equal av- 
erage density. Weinberg's argument, based on photon 
flux conservation, is effectively the basis for the stan- 
dard perturbative approach - i.e. that when averaged 
over distances and angles, the divergence from vacuum 
or underdense regions is compensated by the convergence 
due to clumping, so that the average luminosity distance 
is the same as the luminosity distance in the FL back- 
ground (see e.g. [SS]). Weinberg's argument has been 
disputed |59||60]. Later work (e.g. [6T1|62]) has not pro- 
duced a definitive answer to the question, in particular 
for the case of the very narrow beams involved in SNIa 
observations. 

In order to properly describe a thin geodesic bundle, 
we need to have a good description of the matter distri- 
bution on the scales of the extension of the bundle, and 
determine how the effect of the inhomogeneities average 
during the propagation of the bundle, with two main is- 
sues in mind: (1) determining the typical amplitude of 
the effect and (2) understanding why the description by 
a smooth universe seems to provide a good description 
and determine its validity. 
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A. Angular distance 



For any quantity X evaluated along the ray bundle 



The Jacobi matrix can be diagonalized by rotations: 

= ( (24) 

where the shape parameters D± are nonzero almost ev- 
erywhere. Their absolute values give the semi-axes of the 
(elliptic) cross-section of the bundle. Once D± are fixed, 
the angles ai^2 are unique at all points where the bundle 
is non-circular. 

The area distance or angular diameter distance is then 
defined as'^ 



Da{v) = Vdet\V{v)\ = ^\D+{v)D^{v)\. 



(25) 



For a bundle converging at the observer, Da relates the 
cross-sectional area A at the source to the opening solid 
angle at the observer. It depends on the 4-velocity of 
the observer, but not of the source. From (23), the null 
convergence is 



1 d 



A, 



(26) 



and (21 1 becomes 

dv"^ 



= - (lap - $oo) ^A. 



(27) 



For Tf,^k^'k'' > we have $00 < by (|15|), so that \a\ 
$00 > 0. Thus 



d^DA 
dv^ 



< 0, 



(28) 



in any cosmological model, as long as the null energy con- 
dition holds, irrespective of the value of the cosmological 
constant. 

In order to compare to observations, we need the rela- 
tion between v and z. We have dz/dv — d(u^kf^)/dv ~ 
k-'X/^iuf'k^) = fc^rV^u^. Now 



u^A,, (29) 



where O is the expansion, df^i, is the shear, lo^^^ is the vor- 
ticity and is the acceleration. In a universe containing 
CDM and baryons with four- velocity m'', and A (where 
radiation is dynamically negligible), we have A^ = 0. 
By Eqs. ^ and ([lO|, we obtain [S] 



dv 



= (1 + ^)' 



A^e^ 



(30) 



^ Note that the terminology 'angular diameter distance' has two 
interpretations: as an area angular diameter distance, as used 
here, and as a linear angular diameter distance which are _D+ 
and D- [68) . 



dX \9 TT / i,\ dX 

- = (l + z)^i7,|(.,e'^)-. 



(31) 



where is the observed expansion rate along the line 
of sight P) . 

i^i|(z,e'^) = ie + a^,e''e'^-A^e^ (32) 

The observed expansion rate is made up of an isotropic 
expansion monopole, an acceleration dipole and a shear 
quadrupole. 



The set of equations ( 22 1 , ( 27 ) and ( 30 ) is the basis for 



analyzing the effect of inhomogeneities. There are four 
physical effects induced by inhomogeneities that need to 
be taken into account: 

area distance modifications due to the difference be- 
tween Ricci focussing (when the rays move through 
a uniform medium) and Weyl focussing (due to the 
tidal effects of nearby matter) ; 

redshift adjustment due to the differences between 
the true redshift of a source and its redshift in a 
smoothed out model; 

affine parameter distortions since inhomogeneities 
change the relation v{z) (this is actually where A 
affects observational relations); 

displacement of the light beam since the ray path is 
shifted sideways by inhomogeneities and so experi- 
ences different Weyl and Ricci terms at the same v 
because it is at a different spacetime point. 



B. Links with weak lensing formalism 

The weak lensing amplification matrix A relates the 
direction of observation. 



d?7° 
dv 



to the direction of the source 



so that 



Da{v) ■ 



(33) 



(34) 



(35) 



Here Da is the angular distance in a FL background, 
and — 61- We decompose A into a shear (71,72), a 
convergence k and a rotation w, so that 



A\^ 



1 - K - 7i 72 - w 

72 + CJ 1 — K -I- 71 



(36) 
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The magnification is given by 



(1 



(37) 



where 5, Sq are the surface areas of the image and source 
{S = Sa/detA). Note that while TZat is symmetric by 
construction (see Eq. ( 11 )) and Sab is also symmetric for 



a bundle converging at the observer, it is not necessarily 
the case for T>ab, which actually cannot be generically 
symmetric; see e.g. Eq. (19). 



The amplification matrix ( 35 1 and the deformation ma- 



trix ( 19 ) are both related to the Jacobi matrix, and hence 



they are related by 
av 



A\ — Da 
av 



DaA\S' 



(38) 



This implies (away from caustics, where det^ — 0), 



{A-rM'»y+^b§^=^''>' 

with a prime denoting d/dv and where 



(-4 



_1N a 



1 - K + 7i -72 + UJ 
-72 - w 1 - K - 7i 



The Sachs optical scalars are then given by 

61 = -/i[(l - k)7^ +7ik' + 72^' -a;72] , 

62 = A* [(1 - k)72 +72^' - 71a;' + ^7^] , 



(39) 



(40) 



(41) 

(42) 
(43) 



with the constraint that CJ72 — 721^' = (1 ~ ^)7i ~ Ti'^' 
that arises from V^^fc,^] = 0. These relations are useful 
since the optical scalars are more general (they are de- 
fined for any spacetime geometry) , while the weak lensing 
scalars are widely used in cosmology (but they assume a 
FL background, see Ref. [121 for a general description). 
While generically ui 0, one can see that these equa- 
tions imply that for a FL spacetime only 6 = D'j^/Da 
is non-vanishing at the background level while the shear 
appears only at linear order in perturbation and one has 
(o'l, <5'2) = (^7i7 72) and the rotation appears only et sec- 
ond order in perturbations. 



C. From afHne parameter to redshift dependence 



The evolution equations (22), (27) for the null shear 



and angular distance are in terms of the unobservable 
affine parameter v. We need to convert to the observed 
redshift, using (30). Using (31), we obtain 



= -^(i + z)3i/||e- i(i + z)3fc"v„e 

+ (l-Hz)/c^fc"Vc«A^ + (l + z)^iJ||fc''A^ 
k^ k k V Q (T . 



(44) 



The last term can be evaluated by expanding k" 



cF^vA'^ and 



with (10) and using u^u^Va^u = 
u^e°'\/a^^ = —a^^e°'Vo,u^. It follows that for any quan 
tity X, 



(45) 



— = (1 + zfRl— + (1 + zfQ — , 
dv^ " dz-^ dz 

where 



^/2 



(46) 



This form of Q is completely general, for any space- 
time geometry and energy-momentum tensor, and inde- 
pendent of the field equations. It is convenient to write 
-ffy and Q in terms of covariant multipoles, using a covari- 
ant generalization of a spherical harmonic expansion [64] . 
We expand in terms of the trace-free products e'^^e^\ 
g(Pgj/ga> g^jjjj g(;igi/gag/3) ^ spatial covariant 

derivative V^. Then we obtain |65) : 
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2 4 

-QA^, + -A''a^,^ 



2^ 4 
2e<'^e'^e"U^a,„ + e^^e^'e^e'^V^.a^^, (47) 



and 



Q = 



^(P + 3p) - - A + -e^ + a^.a'^ 
i^^.c.'^- + A^A^ - ^V^A^ 



4 17 

A„ - ^QA.. - -^A^a 



A ^ 



3 ■ 

-I- w^^Wq^ + 2a;QpCT^" — 2Vpv4^ 



(48) 



where we also used the covariant evolution and constraint 
equations of GR (see [66] )• Here £^^,y = C^ai^pu^u^ is the 
electric part of the Weyl tensor (generalizing the Newto- 
nian tidal tensor). These expressions show clearly the 
covariant monopole and higher multipoles; for example. 
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the octupole of Q is V(^<Ti,q,) — ^(^(Ti^q). Note that the 
monopole of iJy has contributions from the shear even 
though the monopole of Hu does not. 



Finally we can rewrite the evolution equation (27) for 
the angular distance in terms of redshift: 



dPA 
dz 



DA.m 



This is a completely general and nonlinear equation, valid 
in any spacetime, with any matter content, where Hu and 



Q are given by (47) and (48). The null shear terms arc 



given by the remaining Sachs equation (22); in terms of 
redshift, this is 



The form of ( 50 1 is unchanged 



IV. 



MODELS BASED ON DIFFERENT 
APPROXIMATIONS 



We briefly review the standard FL approach and the 
DR approximation, and then we propose and investigate 
modifications of the DR model. (For other related re- 
views, see also [67H70].) The set of equations (21), (22) 
and (30) - equivalently (|49|) and (50) - is completely 



general and does not depend on the choice of a particu- 
lar spacetime geometry. We show here how they lead to 
different expressions for the angular distance as a func- 
tion of redshift, depending on the assumptions on the 
distribution of the matter. 



Hu 



dgg 

dz 



1 dDA 



Da dz 

= - {E^, - e^^pe'^H/) N^" , (50) 
Nr ^ {ntn-( ~ n^^nl, n'^.n^ + n^^n'i) , (51) 

where H^i, = \e^a0C"^ uku'^ is the magnetic part of the 
Weyl tensor - which has no Newtonian analogue - and 
£^Lua is the spatial alternating tensor. 



Equations (49) and (501 form a closed system that 
determines Da and a in terms of z and . In par- 
ticular, we see what is required to determine Da in a 
lumpy universe: the total energy-momentum tensor (i.e. 
p,p,q^,'i:^v), the kinematics of the fundamental four- 
velocity (i.e. O, cr^^, oj^y, j4^), the magnetic and electric 
part of the Weyl tensor, H^^ and E^^. 

In a universe with dust matter (CDM and baryons, 
sharing the same four- velocity) , with dark energy in the 
form of A and where we can neglect radiation (i.e. at late 
times), we have 



P 



A,. 



9m 



0. 



(52) 



From now on wc will make this assumption, together with 



0. Then Hy^ and Q simplify to 



Q 



-e^ H — _ 
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fJ.V 



(53) 



AttG 



1 



-P 



-A 



1 



(54) 



and the angular distance equation ( 49 ) becomes 
,d^D, 



(1 + zY H 



+ [1 + z)Q 



dz2 

inGp 



dz 
Da. 



(55) 



A. Smooth FL model 

If we assume the matter is smoothly distributed, then 
the universe can be described by a FL geometry, 

ds' = a'irj) [-dr;2 + dx' + /i(x)dr!'] , (56) 

H{z) = i/o Vf^mo(l + zY + ^AO + ^K0{1 + ^)2(57) 



where fxix) — s\n-{-\/ Kx) / V K is the comoving angular 
distance. The Weyl tensor vanishes, so that — 0, and 
a 

and H\\ 



by (|23|), consistent with (|22|.^lso, Tl^ = $oo(5a, 



( 49 ) reduces to 



H by (31). Then using (15), it follows that 



d^DA 

dz2 



d\nH 
dz 



-2f^mO;^(l 



dDA 
1 + z J dz 

z)Da. 



(58) 



It is important to realize that H{z) in this equation ap- 
pears from the change of variable from v to z. 

This equation also follows directly from the Jacobi ma- 
trix equation ( 17 ) , which is easily solved after a conformal 
transformation, 



dv^'^' 



(59) 



where v is the affine parameter in the conformal space- 
time of the static metric dP = —dif + f'^{x)d^^ ■ The 



solution of (59) is = /k(w)<5^. One can choose v ei- 



ther as 77 or X and the angular distance is then given by 
the standard formula 1 7 1,-73 



DA{z) = j^^fK[x{z)] 



(60) 



Along the past lightconc d^ = —dr] and dz/drj = 
—aoH{z), so that 



aoHoxiz) 



dz' 



H{z')/ Ho- 



rn 
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Since dx = dz/H, (591 recovers (58), after using H 



— (1 + z)HdH/dz and the Einstein equation for H. We 
can either solve (591 or (58) but the first is more direct. 



Also, the derivation of (58 1 required the relation v(z), or 
equivalently v{z). 



For a FL universe Af, 



5t (i.e. 



71 



72 = 0) 



and then 5^ = ifK/fK)^^, so that only 9 = f^^djK/dv 
is non-vanishing, which is a consequence of the spatial 
homogeneity and isotropy. After integration of (22), 
y/A = fx in the conformal spacetime, so that again we 
recover the same expression for the angular distance. 



B. Perturbed FL model 

The simplest way to account for inhomogeneous mat- 
ter is via perturbation theory. At first order for scalar 
perturbations, 



ds' 



a^iv) [-(1 + 2$)d772 + (1 - 2*)7,jda;Ma;Jl , (62) 



in Newtonian gauge, where <i> and 4' are the Bardeen 
potentials. The angular distance is 



Da^Da{1 + Sa), 



(63) 



and the distance duality relation implies that the lumi- 
nosity distance is Dl = (1 + z)'^Da{1 + Sa)- Thus 



SLiz,n) = SA{z,n) + 2 



5z{z, n) 
(1+^ 



2 ' 



(64) 



where n is the direction of observation. The second term 
encodes fluctuations of photon energy due to the local 
gravitational potentials as well as Doppler effects, and 
is similar to the Sachs- Wolfe effect on the CMB. It was 
investigated in [71] and also estimated in [75] in another 
context, but neglected in [TSl [TBI 1751475] . which assumed 
Sl = Sa- 

As long as the bundle remains in the weak lensing 
regime, 



1 + 2k, 



(65) 



so that (5l(z, n) = — ^^(z, n). This assumes that the inho- 
mogeneities can be described by density fluctuations of a 
homogeneous field. Then is a stochastic field of zero 
mean, so that (/i) — 1 (in terms of ensemble average) and 
thus {DA{z,n)) — Da{z). In such a description, the FL 
angular distance is the mean distance that a collection 
of observers will determine. There may indeed be a bias 
from this prediction arising from our actual position in 
the universe. This is a cosmic variance problem. 



1. Derivation from the Jacobi map equation 



The Jacobi equation ( 59 ) reduces to 
d2 



(1) 



T) 

dw2 '^b 



KV 



(1) 

ah 



(66) 



where V^b = vf^ + 2?i',\ and = fK{i)5ab is the 
background Jacobi matrix derived above. This equation 
has the integral solution 



T^abiy)= I fK{i')fK(i-i')ni\\i')di', (67) 



so that the amplification matrix is given by 

N r fK{v')fK{v-v') ,1) , , 



fK{v) 



where 2TZ^^fl — Sg^^^ajsk^'k^n'^n'^ . Since we are inter- 
ested in modes smaller than the Hubble scale, we as- 
sume that the spatial curvature does not influence the 
perturbations and neglect it in the computation of Ti-l^J 
while we keep it in the geometrical factors. We find 
7^^^J,^ = -dadbi^ *) = -29^56$, where = $ since 
we can neglect anisotropic stress at late times and on 
sub-Hubble scales. Then the amplification matrix is 
Aab Sab - dadbip{n, x) with 

r - ^'h [fK{x')n,x'W- (69) 

Jo Jk{x) 

We conclude that (5a = with 

3 

5 A = -K(n,x) = -2^o^mO X 

fK{x')fK{x ~ X') S[fKix')n, x'] , , 

Mx) a(x') ^' ^'"^ 



I 



where the Poisson equation has been used to replace V^^ 
with S. This gives the fluctuation of the angular distance 
in the direction n, taking into account propagation in a 
perturbed spacetime. 



2. Derivation from the Sachs equations 
The same result can, in principle, be derived from the 



first order equations are: 



Sachs equations (21), (22) and (30). The background and 



D 



(0)" 



^00 ' 



.(0)' _ 



r,(i)" _ d,(0)r,(i) , n(0)s(i) 
- ^00 ^A + ^A ^00 ' 



1 + 



lc)(l) , (1) A" 



' 



Sl = Sa + 2- 



.(1) 



l + z(o)' 
where primes are d/dv. 



(71) 
(72) 
(73) 
(74) 

(75) 
(76) 
(77) 
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By (30) with H\\ = H, (72 1 reduces to the definition 
H = aja, and then (71 1 caiTDe integrated to give (60). 



Then ( 73 ) is similar to ( 66 ) , with the source term given 
by 



(1) 

00 



(0) 



1 + z 



(0) 



,(0) 



.(1) 



riO) 



(78) 



Since 9^°'^ is known, (76) can be integrated once the 



source (which depends on gradients of the gravitational 
potentials) is known. 

These two derivations have to give the same answer, 
which they actually do if the effect of perturbations is 
not neglected in any of the equations, and in particular 
in the equation for v{z). 



C. Dyer— Roeder approximation 



of large-scale structure [SSI [HS] . A novel use of the DR 
approach was proposed in [TTl [12] to correct SNIa obser- 
vations for the matter distribution along the line of sight, 
which has important implications for parameter estima- 
tion [86J. 

On large scales (Mpc) we expect a distribution of 3D 
compensated voids giving partially ID compensated mat- 
ter distributions along the line of sight [5S]. The lensing 
effects will not be the same as a smooth distribution of 
matter. On smaller scales (2 kpc to 1 A.U.) we expect 
to mainly move through voids, contaminated by rem- 
nant baryonic gas and non-baryonic dark matter, plus 
the nearly-smoothly distributed photons and neutrinos. 
This suggests that the effect can be significant if matter 
is clustered on small scales with most of the light beams 
used in SNIa observations preferentially moving through 
voids. 



It is important to stress that the perturbative descrip- 
tion still assumes that the distribution of matter is con- 
tinuous (i.e. it assumes that the fluid approximation 
holds on the scale of interest), so that, as long as one is 
in the weak lensing regime, the whole effect arises from 
Ricci focusing with the density of matter equal to the 
average cosmic density (because the effect of the shear 
appears only at second order). Zel'dovich [2] pointed 
out that light is actually more likely to propagate in un- 
derdense regions so that an overall demagnification was 
expected. This idea was followed up by [JHSj. The ap- 
proach came to be named after the later work of Dyer 
and Roeder [SiHSB]. 

The main assumptions of the DR approximation are: 
(1) the Sachs equation ([22]) holds; (2) the relation v{z) 
is the FL one, (30) with i/y = H; (3) the null shear 
a vanishes, as in a FL universe; (4) <I>oo is replaced by 
a{z)^00y where q:{z) represents the fraction of (the mean) 
matter intercepted by the geodesic bundle. In summary, 
the DR model assumes that the bundle is propagating 
in a FL universe but that the Ricci focusing is reduced 
(to reproduce the fact that the beam propagates mostly 
in vacuum or underdense regions) and the Weyl focusing 
remains zero. This implies that the DR equation is 



dz2 



dlni? 
dz 



dPA 
l + zj dz 

{l + z)a{z)DA. 



(79) 



This attempts to model the global effect of inhomogeneity 
in a "mean way" , while still assuming that the universe 
is isotropic and homogeneous. The consistency of the 
DR approximations has been questioned by Ehlers and 
Schneider [Slj and others (e.g. [BTHSi?] ). independently 
of Weinberg's photon flux conservation argument [21] . 

The smoothness parameter a was initially assumed to 
be constant [54H56] . Later it was refined to take into ac- 
coimt its redshift dependence due to the growth of struc- 
ture [521 - IM] and then related to the statistical properties 



D. Modifying the DR approximation 

The previous discussions make it explicit that the DR 
approximation is not a satisfactory model for the effects 
of clumping on ray bundles. Consider the Da{z) equa- 
tion, in the absence of pressure and null shear: 



d 



{l + zfH\\{z) — DA{z) 
= -At:Gp{z)Da{z), 



(80) 



where H\t^{z) = 0/3 -I- a^,ye^e" . The DR approxima- 
tion assumes that we can model the encounters of pho- 
tons with inhomogeneous matter, and not a homogeneous 
background spacetime, by using p{z) as the "true" den- 
sity along a ray - while leaving the rest of the equation 
as in the smooth background. Another way to think of 
this is that the relationship between the affine param- 
eter and redshift is held smooth, and inhomogeneities 
are assumed not to affect the v{z) relation significantly. 
However, photons only experience the local curvature, 
shear and expansion and the average FL-behaviour must 
somehow emerge from integration along the line of sight. 

As discussed above, even in the perturbed FL case the 
relation v{z) fluctuates, and so the DR relation cannot 
be relied upon as a useful approximation even in that sit- 
uation. In particular, a must depend on the line of sight 
since each bundle experiences a different matter profile. 
In a more general sense, the DR approximation does not 
account for changes in the local expansion rate due to 
clumping [SS], and does not capture the essence of weak 
lensing unless a{z) is tuned to a specific form, with no 
apparent physical motivation |85j . 

We do not aim to provide a detailed analysis of the DR 
model here. Rather we offer some possible alternatives in 
order to estimate how significant the effect of clumping 
could be, as well as to show how difficult it is to model 
in a simple but reliable way. 
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1. Modified DR 

In a universe with irrotational dust and arbitrary in- 
homogeneity, the generahzed Friedmann equation is [66) 



1 



-e^ = 87rGp + A 



1 p 



1 



(81) 



where TZ is the Ricci curvature scalar of the 3-surfaces or- 
thogonal to the matter four- velocity. By holding 9 fixed 
to the FL background value 3iJ, the DR approximation 
effectively assumes that the variations of p on any null 
geodesic are compensated by corresponding fluctuations 
in the shear and curvature, which seems unphysical. 

We expect a photon in the real universe to react to the 
nonlocal part of the gravitational field created by dark 
matter halos through local curvature fluctuations in ad- 
dition to the dynamics of the matter in the intervening 
space. A reasonable alternative to DR, then, is to first 
write out the Da{z) equation in a general FL model, 
using the Friedmann equation to evaluate dH/dz. Sub- 
stituting for H{z) using (57) everywhere, we see that 
appears in several places. Then, replacing apm 
gives a plausible alternative to the usual DR approxima- 
tion: 



d'DA 

dz2 



where 



[3a(z)r!„,o(l + z) + nKo] (82) 
^1^^ -2^.0^(1 + .)a(.)i^. 



H{zf - [a{z)n,r,o{l + z)^ + + ^Koil + z)^] 

(83) 

This modified DR equation attempts to take into account 
some aspects of the change in expansion expected from 
an inhomogeneous matter distribution. There are clearly 
a variety of ways to do this (e.g., we have ignored a^z 
terms which could be important), but we have chosen 
just one. See [57] for an alternative approach. 



2. Shell approximation 

Consider a single line of sight, smoothed over some 
scale A. The density profile along this line of sight takes 
some form px{z). If we neglect the angular part of the 
shear cr^^, then we can think of the beam as passing 
through shells of differing density. There exists a spher- 
ically symmetric Lemaitre-Tolman-Bondi (LTB) model 
with non-zero A that has the same density profile p\{z) 
along the past lightcone from the centre. The Da{z) 
relation in the LTB model viewed from the centre will 
approximate the Da{z) relation along the line of sight 
we are trying to model. Each line of sight would have a 
different associated LTB model (a mosaic of cones around 



us, in the language of [H]). The utility of this approxima- 
tion lies in the fact that we can specify a density profile 
on a surface of constant time, and use the exact LTB 
solution to evolve the density backwards onto the past 
lightcone. We can then calculate Da{z) exactly for that 
line of sight. Most importantly, this will account for the 
variable expansion rate along the direction of propaga- 
tion which also takes into account the radial component 
of the shear. 

In the LTB model, the angular Hubble rate is given by 
an effective Friedmann equation [88] 



= nn,Q{r) aj^ + nKo{r) + l^Ao(r), (84) 



where the angular scale factor a±{t,r) is normalized to 
unity today, the i7's have an arbitrary radial degree of 
freedom in them, and H±o(r) is calculated once the age 
is set (or the Hubble rate at the centre is chosen) . If fiuio 
is chosen as a constant we have an FL model. To model 
a radial line of sight we can choose the density profile 
today as po{r) = [1 + S{r)]px{0), and then 



1 



HUr) 



dr'r' pa{r') . 



(85) 



Then (84 1 evolves the density back onto the past light- 



cone, using 

di 
d^ 



1 



{l + z)Hn 



where the radial 

[dtdr{a±r)]/[dr{a±r)] 



dr 



flKoHl^r'- 



(1 + z)dtdr{a±_r) 



(86) 



Hubble rate is i/y (i, r) = 
The area distance is then 



DA{z)^a^{t{z),r{z))r{z). 



3. Numerical investigation 



(87) 



We can compare these different approximations numer- 
ically. First consider the case of a single density fluctua- 
tion. Figure|4]shows the results of looking through a large 
void and large overdensity, 500 Mpc away, modelled with 
a Gaussian deviation from a = 1 of width ^ 100 Mpc. An 
underdensity causes an increase in the distance modulus 
at redshifts beyond itself, in both the DR and modified 
DR cases; the opposite happens for an overdensity. 

The DR and modified DR are qualitatively similar, 
while the shell approximation is very different. Accord- 
ing to the (modified) DR approximations, we should ex- 
pect SNIa to appear dimmer when located behind an 
underdensity as compared to an overdensity. By con- 
trast, the shell approximation gives the opposite effect 
with a much larger amplitude: SNIa located behind a 
void appear brighter than in the fiducial cosmology; lo- 
cated behind an overdensity, they appear dimmer. The 
reason is as follows: although a void results in a nega- 
tively curved region (which would imply diverging light 
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shell 

FL background 




FIG. 4: The effect of a single void (left) or overdensity (right) 
on the distance modulus fi = 5 logj^Q(_D_L/10 pc) according to 
three different approximations, using as a 'background' a flat 
LCDM model with Qm = 0.25, h = 0.7. 



rays and larger distances), this is accompanied by an in- 
crease in the expansion rate in that region, which actually 
has a much stronger effect on distances (compare |69j). 
In FL, increasing the expansion rate and decreasing the 
density while keeping the age fixed results in a model with 
smaller distances, and this is exactly what happens here. 
For an overdensity, the reverse applies. However, note 
that while the shell approximation captures the mean 
expansion rate down the line of sight nicely, it may not 
capture the radial expansion rate correctly. E.g., looking 
through a spherical void vs a shell with the same density 
profile have different shear along the line of sight, making 
our approximation somewhat exaggerated in such a case. 

Now consider the case where the density along the line 
of sight is reduced by a fixed percentage below the back- 
ground value. Figure [3] shows that the main contribution 
of smoothing a simulation over smaller beam sizes is to 
reduce the mean value of a, since particles of significant 
size are rarely encountered. In the DR and modified DR 
approximations this amounts to fixing a = const, and 
Monte Carlo-ing over the PDF. In the shell approxima- 
tion we have S = a — 1 (since p is constant, this is re- 
ally just an FL model with adjusted parameters). In 
Fig. [5] we show the distance modulus for two fiducial 
backgrounds: one a standard LCDM model, the other 



possible change to the distance modulus 

0.6] ^ approximation: 

DR 

niDR 

sheU 




FL 

background 



0.1 



1.0 



FIG. 5: We show the effect of differing mean values 
of a based on different approximations described in the 
text. The best known of these is the Dyer-Roeder (DR) 
which simply reduces the energy density in the light prop- 
agation equations while keeping the background expan- 
sion rate. Attempts to model this more accurately by 
accounting for the varying expansion rate encountered 
along the beam yield very different contributions depend- 
ing on how this is modelled. One such modified DR 
(mDR) gives dark energy-like behaviour in the distance 
modulus even for a model with no cosmological constant 
(all shown compared to an empty model). 



a curved CDM model with A = 0. We see that all the 
approximations give a systematic dimming effect to vary- 
ing degrees. (There is actually an overall brightening in 
the shell case due to the change in h, but we have sub- 
tracted this off because SNIa observations are only sen- 
sitive to relative magnitudes, not absolute ones, which is 
marginalised over together with h.) While the DR and 
shell approximations are rather similar, giving changes to 
the distance modulus of at most 0.1 at z ~ 1, we see that 
the modified DR approximation gives a much stronger 
effect, of several times that. 

It is striking that in the case A = for a < 0.3 we see 
the modified DR distance modulus mimicking the be- 
haviour of a dark energy component. More specifically, 
both the DR and modified DR approximations can mimic 
a LCDM model, though the DR approximation requires 
a drastic change to a at low redshift to do so. For the 
modified DR approximation this is not so - see Fig. (|6| - 
and it is well known that an LTB distance modulus can 
mimic any FL one. 

It is clear from Fig. |4] that the DR approximation and 
plausible variations of it can give very different results - 
so that the basis of the DR approximation itself is sus- 
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FIG. 6: The a{z) required to give a LCDM Da{z) curve, for 
a variety of ilmo (with curvature making up the rest of the en- 
ergy component). The DR approximation requires extremely 
negative values of a to mimic dark energy, which might sug- 
gest that the effect of modelling narrow beams may not be 
important for dark energy, and certainly could not be the 
underlying cause of it. On the other hand, a simple modifica- 
tion to DR yields an approximation that much more readily 
produces dark energy-like effects. This suggests instead that 
modelling narrow beams properly, and getting the DR ap- 
proximation right, could be vital for determining the nature 
of dark energy. 



tions (e.g. [55]) is very useful for statistical analysis and 
predictions for lensing observations. However, the dy- 
namical range of scales of matter inhomogeneities that 
the simulations can reproduce is very limited. Ray- 
tracing within N-body simulations is usually done by pro- 
jecting all matter onto equally spaced 'lens planes' which 
are separated by about 100 Mpc. Such an approach is not 
able to clarify the effect that halos below the simulations 
mass resolution have on light bundles from SNIa. 

Exact solutions that are more general than the Swiss 
cheese models are usually less realistic, given the highly 
nonlinear nature of GR. A class of Szekeres models, in 
which inhomogeneities may be modelled as nonlinear per- 
turbations on an FL background, has been used to inves- 
tigate light propagation by [31]. Despite the idealized 
nature of clumping, the results show that for inhomo- 
geneities of large spatial size, parameter estimation could 
be seriously affected (see also I100|). A stronger conclu- 
sion follows from analyzing light rays in a universe with 
regularly spaced point masses separated by vacuum |101j . 
The average dynamics is close to LCDM, but the op- 
tics behaves very differently. Unlike N-body simulations, 
this model is self-consistent (i.e. the particles generate 
their own spacetime geometry). However the modelling 
of matter is necessarily over-simplified. 

V. RICCI AND WEYL FOCUSING 



pect. 

E. Other approaches 

Other approaches to the problem of light propagation 
in a clumpy universe are based on exact nonlinear solu- 
tions or numerical methods or both. 

In "Swiss cheese" models, an FL universe contains 
one or more spherical inhomogeneous regions, with the 
same average density as the FL model, which could be 
Schwarzschild vacuoles, or more realistically, LTB balls. 
These have been used extensively to model the effect of 
inhomogeneities on light rays (e.g. [SSH^ I. The results 
depend on the position and nature of the inhomogeneous 
regions, as well as on the method for randomization of 
light rays. A careful analysis that approximates obser- 
vations in all directions and over a range of distances 
[97] concludes that there are only small corrections to 
the standard FL results. This is not surprising, since 
the average density along typical lines of sight is very 
close to the FL density. For the case of SNIa beams, 
which may preferentially sample underdense regions, the 
results could be different. However, an inherent limita- 
tion of this approach is that, by construction, and given 
the highly symmetric nature of the inhomogeneities, the 
clumps do not affect the expansion rate of the universe. 

The method of ray-tracing through N-body simula- 



The DR approximation neglects the effect of point 
sources and the Weyl focusing they produce, in particu- 
lar in the strong lensing regime where it is not negligible. 
This issue was addressed in [21], by considering the ef- 
fect of the DR equation (for almost all directions, where 
Weyl effects can be neglected) and the effect of strong 
lensing (for the relatively few directions where lightrays 
pass close to matter, with strong lensing occurring and 
leading to multiple images). These combined effects are 
shown to lead to the usual FL relations when averaged 
over the whole sky, the decreased flux in most directions 
being compensated by higher flux in a few directions. 

Keeping in mind that strong lensing effects are negli- 
gible in most directions, and in particular for most SNIa 
(unless we observe a galaxy or a cluster on the line of 
sight), we can try to go a step further than the DR ap- 
proximation by relating, at least heuristically, the Weyl 
focusing to an effective Ricci focusing. The approxi- 
mation that the SNIa bundles remain in the weak-field 
regime can be supported as follows. If matter is mod- 
elled as a gas of particles of mass m and proper radius 
r*, with mean number density n, then the mean energy 
density is p = mn. The probability that a line of sight 
intersects such a matter particle within the redshift band 
z — > z-|-dz is proportional to the surface area of the parti- 
cles, to their density, and to the distance light propagates 
in this redshift interval, i.e. 
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The average number of particle intersections before the 
redshift z is the optical depth, 



r(z) 



dz'. 



(89) 



We assume particle number conservation, n(2;) = '^o(l + 
z)^, with no ^ 0.005/i'^Mpc^'^ (number density of ha- 
los above 10^^ h~^MQ which comprise about 50% of the 
mass in the universe) and 20/i~^kpc for the central, 
high density, region of the halos where the galaxies are 
assumed to reside. This implies that at z = 1, r ~ 0.023 
and T = 0.032, which means that only 2.3% and 3.2% of 
the lines of sight intersect a galaxy before z — 1, respec- 
tively for an Einstein-de Sitter and a flat ACDM model 
with r2,„ = 0.3. This is a rough estimate - A and in- 
homogeneous distribution of the luminous matter would 
tend to lower it. 

In order to describe the transition from Wcyl to Ricci 
focusing, we first recall the lensing effect of a single mass, 
whose gravitational potential is 



-Gm{b^ 



-1/2 



(90) 



where b is the impact parameter and u the distance along 
the line of sight. The deflection angle is thus 



4Gm 



(91) 



The critical density and the Einstein angular radius are 
respectively defined by Scrit — Dqs / {^t^GDolDls) and 
= ^GmDLs I {DolDqs), so that the lens equation 
takes the form 



e = e^ + a.{e), 



(92) 



where D^s, Dql and Dqs are the angular distances be- 
tween the lens plane and the source, the observer and 
the lens plane, and the observer and the source. The 
angular position of the image on the lens plane is then 
given hy — b/DoL- ct is given by a = 60^/0'^, since 
E(6») = m(5(2)(5) = m(5(2) (6/)/d2^^ for a point mass. The 
amplification matrix is then obtained as A^^ = 89^ /d9^, 
so that Aab = Sab - Sgadb = Sab ~ de^debijj with 



9|ln( 



(93) 



We conclude that the convergence and the shear are 
given by 

^ = n0l5^'\e), (71,72) = §(^?2'-^?, 20102), (94) 

which is more easily written in terms of the polar angle 
on the screen (0i = 9 cos (p, 02 = 9 sin (p) as 



(71,72) 



02 



(cos2(^, sm2(p). 



(95) 



For a single mass and a line of sight that does not 
intersect it, k = and the magnification is fi = [1 — 



|7|2]~i ^1-1- 20|0~4^ neglecting the image inside the 
Einstein radius. Now consider a shell of thickness dz 
such that the density of particles is n{z) and such that 
dx ^ Dls, DoL, Dos- For a typical lightray the total 
amplification matrix will have a shear given by 



[71 (^), 72(2:)] 



2 v;^ rcos2^ 
^e2^[ 02 ' 



sin 2ip 



02 



(96) 



where the sum is over all particles of the shell. If the 
particles are distributed homogeneously, this implies that 
they are isotropically distributed around the line of sight, 
so that we expect 71(2:), 72(2:) ~ 0. 

Intuitively, this is understood by the fact that each 
point mass induces an ellipticity in a different direction 
and they should average to reduce the total ellipticity. 
The remaining effect of all the Weyl distortions is thus 
an effective convergence that can be determined from the 
magnification k{z) — 02^0^4 Estimating the sum by 
assuming we have a uniform distribution and that the 
typical smallest distance is ~ n^^^"^ /Dql, we get that 



4 Dqj^ti, and thus 



Kiz) - fMxiz)Hz) \ AGm 



yK[x{z)]fK[xizs)] 



(97) 



where fxix) is the unfilled angular distance, i.e. using H 
instead oi H in (61 ). We can compare to the effect that a 
homogeneous distribution of density nes{z) = a{z)n{z) 
would generate through its Ricci focusing to get 



a{z) ~ AGm 



fK[x{zs) - x{z)] 
fK[x{zs] 



H{z). 



(98) 



This is the effective DR parameter for the averaged Ricci 
convergence. It stills depend explicitly on m because this 
is a second order effect, hence scaling as 0g, while it is 
first order in a homogeneous medium. Such an estimate 
is indeed very crude but it confirms the statements of 
f24i |57] and the results of the numerical simulations of 

im. 



VI. DISCUSSION AND CONCLUSIONS 

The effect of the inhomogeneity of the matter distri- 
bution induces a dispersion of the magnification of SNIa 
and thus of the luminosity distance. We have argued that 
this effect has not been properly modelled for SNIa since 
the beam is very narrow and far below the scales resolved 
in any numerical simulation. 

For the first time, we have attempted to quantify the 
probability distribution for narrow beams, using a com- 
bination of N-body simulations and a PS approach. For 
a narrow beam of fixed length, the PDF is non-Gaussian, 
peaked at densities below the cosmic mean, with a power- 
law tail, whose power depends on the diameter of the 
beam, describing the relatively few lines of sight which 
have an over-dense mean. These PDFs contrast sharply 
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with distributions based on using cubes of the same vol- 
ume. These estimations are based on current N-body 
simulations which do not have the resolution to probe 
beams with a diameter < 100/i~"'^kpc. Nevertheless, the 
trend is clear: narrow beams typically experience a lower 
than average density, and do not sample the cosmic mean 
density until their length approaches the Hubble scale. 
Based on our results, we estimate that significantly more 
than 75% of beams experience less than the mean density. 
From a theoretical point of view, the effect must be 



described by the set of equations (21 ), (22) and (30) that 



describe the distortion and magnification of any light 
bundle, whatever the spacetime geometry. The explicit 
covariant form of these equations is given by (50)~(55|. 



The main problem is that the solution of this set of equa- 
tions requires a description of the distribution of matter 
on the scale of the beam size, i.e. on scales much smaller 
than those of our current understanding. 

This dispersion has been modelled in some regimes but 
we argue that the distribution of matter on the scales rel- 
evant for the description of a SNIa bundle is not under- 
stood yet. On such small scales, the statistical dispersion 
comes with a bias that has two origins: (1) the fact that 
the nonlinear terms in the expression for the magnifica- 
tion can not be neglected a priori; (2) an observational 
selection effect due to the fact that most SNIa are ob- 
served in directions where they are not overshadowed by 
a galaxy. The bundles included in SNIa analysis are thus 
more likely to probe under-dense regions. 

Why is the Hubble diagram so compatible with that of 
a Friedmann-Lemaitre universe? In particular, the typi- 
cal transverse size of the bundle is smaller than the typi- 
cal mean distance between the smallest bound structures 
in the currently favoured CDM paradigm. Why does 
the fluid approximation used to interpret the data hold? 
This question is two-sided. It questions the robustness of 
the interpretation of the cosmological data but also of- 
fers a way to constrain the distribution of matter on small 
scales. We also gave a heuristic argument concerning the 
Ricci-Weyl focusing issue, leading to a prediction of an 
effective DR factor a{z) once the fractions of clustered 
and smooth matter are known. 

Another description is provided by the Dyer-Roeder 
equation. It has however some simplifying hypotheses 
that neglect the effects of changes in v{z) and the local 
expansion due to clumping. We suggested two plausible 
modifications to the DR approach, and showed that the 
3 models produce very different results - thus undermin- 
ing confidence in the DR approximation. In particular, 
it is not clear whether under-densities lead to demagnifi- 
cation (due to negative curvature) or magnification (due 
to the increase in the expansion rate) . Our shell approxi- 
mation clearly points to the opposite effect calculated via 
normal lensing or the DR approximation: a SN located 
just behind an under-density should appear brighter than 
it would in the fiducial cosmology. In fact, we estimate 
that a SN the far side of a 100 Mpc void could be 0.1 mag 
brighter than it would be with not void present. Though 



likely an over-estimation, this could have important im- 
plications for parameter estimation from SNIa. Quanti- 
fying this properly is an important open problem. In fact, 
it is striking to note from our investigations of N-body 
simulations that we should expect a to vary as a function 
of radius from us from significantly below unity locally, 
approaching unity only on Hubble scales. Within our 
shell approximation, this would give exactly the kind of 
model used in Hubble scale 'void models', which require 
no dark energy, but without any kind of anti-Copernican 
fine tuning involved |103| . Every observer would observe 
such an effect. 

While accurate modelling of such beams may be prob- 
lematic for some time, we can still observationally test to 
see if there are problems and phenomenologically correct 
for them. 

SNIa line of sight Dividing up SNIa samples accord- 
ing to the estimated density along the line of sight 
may reveal a bias. If so, this may indicate that the 
effects we have discussed here must be taken into 
account. 

discordant distances In any exact relativistic model 
the luminosity distance is Dj^ = (1 -I- z)^Dji, where 
Da is the area angular diameter distance. Large 
scale measurements of the area distance will not be 
affected by the problems we have discussed here, 
however, so we can expect a failure at some level of 
this relation when comparing measurements of Da 
from large scale measurements such as the BAO 
and the CMB to measurements of from SNIa. 
On smaller scales, where the area distance is mea- 
sured from radio or quasar sources, there could 
be an effective reciprocity breakdown because such 
sources still have much larger beam sizes than SNIa, 
and so smear the matter distribution to include 
many more over-dense regions. Such a violation 



was found in [102j . where a relative brightening of 
SNIa was found - as we would predict from the 
shell approximation. 

consistency conditions A variety of consistency tests 
have recently been developed as a way of testing 
the standard model (see |103j for a review). It has 
recently been shown that these are strongly sensi- 
tive to changes of the DR form pII4] . This implies 
that they can be used to probe the effects we have 
discussed here. 

Generically, then, distances to the same object will de- 
pend on the scale over which the light from the source 
smears the intervening matter distribution. 

We have found that the old problem of modelling nar- 
row beams remains unsolved. As different interpretations 
of the problem give conflicting yet significant effects, we 
believe this problem needs considerably more attention. 
This is important not only from a theoretical perspec- 
tive, but to ensure precision cosmology delivers correct 
answers as well as precise ones. 
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